#################################
###   Replication Code for    ###
### Studying Malapportionment ###
### Using <alpha>-Divergence  ###
###       Apr. 20, 2018       ###
###      Yuta KAMAHARA        ###
### (Yokohama National Univ.) ###
#################################
# FOR WINDOWS
# Environment: 
# Windows 32bit: Java 32bit
# Windows 64bit: Java 64bit
# FOR MAC OS
# Users can directly download Samuels and Snyder's (2001) data set from the following URL
# ("http://images.cla.umn.edu/samuels-snyder-data.zip")
# and then calculate our measures using Sheet "Lower Chambers" in their data set.

library(xlsx)
library(ggplot2)
library(gridExtra)
library(reshape2)

################
### Figure 2  ##
################

# Lorenz's dominance
#45 degree
fortyfive <- cbind(
  "Accumulated Ratio of Population" = c(0, 0.25, 0.5, 0.75, 1),
  "Accumulated Ratio of Apportionment" = c(0, 0.25, 0.5, 0.75, 1)
)

# Case1
case1 <- cbind(
  "Accumulated Ratio of Population" = c(0, 2000/6000, (2000/6000)*2, ((2000/6000)*2)+(1000/6000), 1),
  "Accumulated Ratio of Apportionment" = c(0, 1/6, (1/6)*2,((1/6)*2)+(2/6), 1)
)

# Case2
case2 <- cbind(
  "Accumulated Ratio of Population" = c(0, 3000/6000, (3000/6000)+(1000/6000), (3000/6000)+(1000/6000)*2, 1),
  "Accumulated Ratio of Apportionment" = c(0, 1/6, (1/6)*2, (1/6)*3, 1)
)

lbllrz <- c("45 Degrees", "Case 1", "Case 4")

pdf("Figure2.pdf")
par(xaxs="i", yaxs="i") #set the origin
plot(fortyfive, type="l", col="grey", lwd=4, panel.first=grid(lty=1))# writing grids
par(new=T)
plot(case1, type="l", lwd=2)
par(new=T)
plot(case2, type="l", lty=2, lwd=2)
legend("topleft", lbllrz, col=c("grey", "black", "black"), lty=c(1,1,2), lwd=c(4,2,2), bg="white")
dev.off()

###############################
###        loading          ###
### Samuels & Snyder (2001) ###
###############################

loc.url <- "http://images.cla.umn.edu/samuels-snyder-data.zip"
td <- tempdir() 
tf <- tempfile(tmpdir=td, fileext=".zip") 
download.file(loc.url, tf) 

fname <- unzip(tf, list=TRUE)$Name[6] 
unzip(tf, files=fname, exdir=td, overwrite=TRUE) 
fpath <- file.path(td, fname)
data <- read.xlsx(fpath, sheetIndex = 1)
unlink(td)
rm(td)
unlink(loc.url)
rm(loc.url)
unlink(tf)
rm(tf)
unlink(fpath)
rm(fpath)

data[339,7]
data[339,7]<-50369
data2 <- data[,-4]  # comma separation has errors because of Hungary; thus we deleted a var. for district name.
write.csv(data2, "SamuelsSnyder.csv", row.names=FALSE)
ss <- read.csv("SamuelsSnyder.csv", header=TRUE)
table(ss$X.SEATS)
table(data$COUNTRY)
table(data$CODE)
# table(ss$COUNTRY)

# generate new ccode and its definition
ss$country.id <- NA
ss$country.id[ss$COUNTRY=="ANDORRA"] <- 1
ss$country.id[ss$COUNTRY=="ANTIGUA AND BARBUDA"] <-	2
ss$country.id[ss$COUNTRY=="ARGENTINA"] <- 3
ss$country.id[ss$COUNTRY=="AUSTRALIA"] <-4
ss$country.id[ss$COUNTRY=="AUSTRIA"] <- 5
ss$country.id[ss$COUNTRY=="BARBADOS"] <- 6
ss$country.id[ss$COUNTRY=="BELIZE"] <- 7
ss$country.id[ss$COUNTRY=="BENIN"] <- 8
ss$country.id[ss$COUNTRY=="BOLIVIA"] <- 9
ss$country.id[ss$COUNTRY=="BRAZIL"] <- 10
ss$country.id[ss$COUNTRY=="BURKINA FASO"] <- 11
ss$country.id[ss$COUNTRY=="CANADA"] <- 12
ss$country.id[ss$COUNTRY=="CHILE"] <- 13
ss$country.id[ss$COUNTRY=="COLOMBIA"] <- 14
ss$country.id[ss$COUNTRY=="COSTA RICA"] <-	15
ss$country.id[ss$COUNTRY=="CYPRUS"] <- 16
ss$country.id[ss$COUNTRY=="CZECH REP."] <-17
ss$country.id[ss$COUNTRY=="DENMARK"] <- 18
ss$country.id[ss$COUNTRY=="DOMINICAN REPUBLIC"] <- 19
ss$country.id[ss$COUNTRY=="ECUADOR"] <- 20
ss$country.id[ss$COUNTRY=="EL SALVADOR"] <- 21
ss$country.id[ss$COUNTRY=="ESTONIA"] <- 22
ss$country.id[ss$COUNTRY=="FINLAND"] <- 23
ss$country.id[ss$COUNTRY=="FRANCE"] <- 24
ss$country.id[ss$COUNTRY=="GAMBIA"] <- 25
ss$country.id[ss$COUNTRY=="GEORGIA"] <- 26
ss$country.id[ss$COUNTRY=="GERMANY"] <-	27
ss$country.id[ss$COUNTRY=="GHANA"] <- 28
ss$country.id[ss$COUNTRY=="GREECE"] <-29
ss$country.id[ss$COUNTRY=="GUATEMALA"] <- 30
ss$country.id[ss$COUNTRY=="HONDURAS"] <- 31
ss$country.id[ss$COUNTRY=="HUNGARY"] <- 32
ss$country.id[ss$COUNTRY=="ICELAND"] <- 33
ss$country.id[ss$COUNTRY=="INDIA"] <- 34
ss$country.id[ss$COUNTRY=="IRELAND"] <- 35
ss$country.id[ss$COUNTRY=="ITALY"] <- 36
ss$country.id[ss$COUNTRY=="JAMAICA"] <- 37
ss$country.id[ss$COUNTRY=="JAPAN"] <- 38
ss$country.id[ss$COUNTRY=="KENYA"] <-	39
ss$country.id[ss$COUNTRY=="KOREA"] <- 40
ss$country.id[ss$COUNTRY=="LATVIA"] <-41
ss$country.id[ss$COUNTRY=="LIECHTENSTEIN"] <- 42
ss$country.id[ss$COUNTRY=="MALAWI"] <- 43
ss$country.id[ss$COUNTRY=="MALI"] <- 44
ss$country.id[ss$COUNTRY=="MALTA"] <- 45
ss$country.id[ss$COUNTRY=="MEXICO"] <- 46
ss$country.id[ss$COUNTRY=="NEW ZEALAND"] <- 47
ss$country.id[ss$COUNTRY=="NICARAGUA"] <- 48
ss$country.id[ss$COUNTRY=="NORWAY"] <- 49
ss$country.id[ss$COUNTRY=="PANAMA"] <- 50
ss$country.id[ss$COUNTRY=="PARAGUAY"] <-	51
ss$country.id[ss$COUNTRY=="POLAND"] <- 52
ss$country.id[ss$COUNTRY=="PORTUGAL"] <-53
ss$country.id[ss$COUNTRY=="ROMANIA"] <- 54
ss$country.id[ss$COUNTRY=="RUSSIA"] <- 55
ss$country.id[ss$COUNTRY=="SENEGAL"] <- 56
ss$country.id[ss$COUNTRY=="SEYCHELLES"] <- 57
ss$country.id[ss$COUNTRY=="SLOVAKIA"] <- 58
ss$country.id[ss$COUNTRY=="SLOVENIA"] <- 59
ss$country.id[ss$COUNTRY=="SOUTH AFRICA"] <- 60
ss$country.id[ss$COUNTRY=="SPAIN"] <- 61
ss$country.id[ss$COUNTRY=="SRI LANKA"] <- 62
ss$country.id[ss$COUNTRY=="ST. LUCIA"] <- 63
ss$country.id[ss$COUNTRY=="SWEDEN"] <-64
ss$country.id[ss$COUNTRY=="SWITZERLAND"] <- 65
ss$country.id[ss$COUNTRY=="TANZANIA"] <- 66
ss$country.id[ss$COUNTRY=="THAILAND"] <- 67
ss$country.id[ss$COUNTRY=="TURKEY"] <- 68
ss$country.id[ss$COUNTRY=="UKRAINE"] <- 69
ss$country.id[ss$COUNTRY=="UNITED KINGDOM"] <- 70
ss$country.id[ss$COUNTRY=="URUGUAY"] <- 71
ss$country.id[ss$COUNTRY=="USA"] <- 72
ss$country.id[ss$COUNTRY=="VENEZUELA"] <- 73
ss$country.id[ss$COUNTRY=="ZAMBIA"] <- 74
ss$country.id[ss$COUNTRY=="ZAMBIA "] <- 74 # unexpected existence of blank space
ss <- ss[!ss$X.SEATS==0,] # Panama had one electoral district with 0 district-magnitude.

ss <- ss[,c(1,2,5:7,11)]
sum_seat<- tapply (ss$X.SEATS, ss$country.id, sum)
sum_seat

sum_pop <- tapply (ss$DISTRICT.POPULATION, ss$country.id, sum)
sum_pop
Sum <- rbind(sum_seat, sum_pop)
Sum <- t(Sum) # swapping rows with columns 
ss_name <- ss[!duplicated(ss$CODE),1:2]
Sum[1,1] # Andorra, col1
Sum[1,2] # Andorra, col2

Sum[74,1] # Zambia, col1
Sum[50,1]
Sum[50,2]
1500451-1496751 # check - Panama: # of electorates was 3700 in Embera

cname <- ss_name[,1] 
rownames(Sum) <- cname
result <- as.data.frame(matrix(NA,nrow=length(unique(ss$country.id)),ncol=8))

for (i in 1:length(unique(ss$country.id))){               # loop for country  
  c <- subset(ss, country.id==i, select=c(X.SEATS,DISTRICT.POPULATION))
  P <- c$DISTRICT.POPULATION / Sum[i,2]
  Q <- c$X.SEATS / Sum[i,1]
  mal <- sum(abs(P-Q))/2
  Adams <- sum(P*(1/((-9)*((-9)-1))) * (((Q/P)^(-9)-1))) # alpha = -9
  Hill <- sum(P*(1/((-1)*((-1)-1))) * (((Q/P)^(-1)-1)))
  Nash <- sum(P* (-log(Q/P))) 
  Bentham <- sum(P * ((Q/P) * log(Q/P)))  
  Webster <- sum(P*(1/(2*(2-1))) * (((Q/P)^2)-1))
  Jefferson <- sum(P*(1/(10*(10-1))) * (((Q/P)^10)-1)) #alpha = 10
  result[i,1] <- mal
  result[i,2] <- Adams
  result[i,3] <- Hill
  result[i,4] <- Nash
  result[i,5] <- Bentham
  result[i,6] <- Webster
  result[i,7] <- Jefferson  # save results
  result[i,8] <- i  
}

colnames(result) <- c("mal", "Adams", "Hill", "Nash", "Bentham", "Webster", "Jefferson", "ID")
rownames(result) <- cname

write.csv(result, "divergence.csv", quote=FALSE, row.names=TRUE)
divergence2 <- result[-50,] # drop Pananma because 0 seat for 1 district
write.csv(divergence2, "divergence2.csv", quote=FALSE, row.names=TRUE)

rankD <- divergence2
rankD <- transform(rankD, r_mal=rank(mal, ties.method="min")) 
rankD <- transform(rankD, r_Adams=rank(Adams, ties.method="min")) 
rankD <- transform(rankD, r_Hill=rank(Hill, ties.method="min")) 
rankD <- transform(rankD, r_Nash=rank(Nash, ties.method="min")) 
rankD <- transform(rankD, r_Bentham=rank(Bentham, ties.method="min")) 
rankD <- transform(rankD, r_Webster=rank(Webster, ties.method="min")) 
rankD <- transform(rankD, r_Jefferson=rank(Jefferson, ties.method="min")) 


cor.test(rankD$Nash, rankD$Bentham, method="spearman")
cor.test(rankD$Nash, rankD$Bentham, method="kendall") #Fig.3
cor.test(rankD$Hill, rankD$Webster, method="spearman")
cor.test(rankD$Hill, rankD$Webster, method="kendall") #Fig.3
cor.test(rankD$Adams, rankD$Jefferson, method="spearman")
cor.test(rankD$Adams, rankD$Jefferson, method="kendall") #Fig.3

write.csv(rankD, "rank.csv" , quote=FALSE, row.names=TRUE)

################
### Figure 3  ##
################

#0.1 rank correlation, Spearman & Kendall, MAL & Nash
cor.test(rankD$mal, rankD$Nash, method="spearman")
mn_s <- expression("Spearman's"~rho == 0.976)
cor.test(rankD$mal, rankD$Nash, method="kendall") 
mn_k <- expression("Kendall's"~tau == "0.890")

cma_mn1 <- ggplot(rankD, aes(x=r_mal, y=r_Nash))
cma_mn1 <- cma_mn1 + geom_point() + theme_bw() + 
  xlab("MAL-based rank")+ 
  ylab(expression(paste("Nash-based rank (", alpha %->% 0, ')')))+
  annotate("text", x=22, y=70, parse = T, label=as.character(mn_s), size=9) +
  annotate("text", x=19, y=65, parse = T, label=as.character(mn_k), size=9) +
  labs(title="(a)")+
  theme(plot.title = element_text(size = rel(2)))+
  theme(axis.title.x = element_text(size=20),axis.title.y = element_text(size=20)) 

cma_mn2 <- ggplot(rankD, aes(x=r_mal, y=r_Nash))
cma_mn2 <- cma_mn2 + geom_point() + theme_bw() + 
  xlab("MAL-based rank")+ 
  ylab(expression(paste("Nash-based rank (", alpha %->% 0, ')')))+
  annotate("text", x=22, y=70, parse = T, label=as.character(mn_s), size=4) +
  annotate("text", x=19, y=65, parse = T, label=as.character(mn_k), size=4) +
  annotate("text", x=46, y=60, label="Canada", size=3) +
  annotate("text", x=34, y=57, label="India", size=3) +
  labs(title="(a)")+theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.title = element_text(size = rel(1)))+
  theme(axis.title.x = element_text(size=12),axis.title.y = element_text(size=12)) 

#0.2 rank correlation, Spearman & Kendall, MAL & Bentham
cor.test(rankD$mal, rankD$Bentham, method="spearman")
mb_s <- expression("Spearman's"~rho == 0.965)
cor.test(rankD$mal, rankD$Bentham, method="kendall") 
mb_k <- expression("Kendall's"~tau == 0.865)

cma_mb1 <- ggplot(rankD, aes(x=r_mal, y=r_Bentham))
cma_mb1 <- cma_mb1 + geom_point() + theme_bw() + 
  xlab("MAL-based rank")+ 
  ylab(expression(paste("Bentham-based rank (", alpha %->% 1, ')')))+
  annotate("text", x=22, y=70, parse = T, label=as.character(mb_s), size=9) +
  annotate("text", x=19, y=65, parse = T, label=as.character(mb_k), size=9) +  labs(title="(d)")+
  theme(plot.title = element_text(size = rel(2)))+
  theme(axis.title.x = element_text(size=20),axis.title.y = element_text(size=20)) 

cma_mb2 <- ggplot(rankD, aes(x=r_mal, y=r_Bentham))
cma_mb2 <- cma_mb2 + geom_point() + theme_bw() + 
  xlab("MAL-based rank")+ 
  ylab(expression(paste("Bentham-based rank (", alpha %->% 1, ')')))+
  annotate("text", x=22, y=70, parse = T, label=as.character(mb_s), size=4) +
  annotate("text", x=19, y=65, parse = T, label=as.character(mb_k), size=4) +
  annotate("text", x=46, y=55, label="Canada", size=3) +
  annotate("text", x=34, y=56, label="India", size=3)  +
  labs(title="(b)")+theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.title = element_text(size = rel(1)))+
  theme(axis.title.x = element_text(size=12),axis.title.y = element_text(size=12)) 


#1 rank correlation, Spearman & Kendall, Nash & Bentham
cor.test(rankD$Nash, rankD$Bentham, method="spearman")
nb_s <- expression("Spearman's"~rho == 0.996)
cor.test(rankD$Nash, rankD$Bentham, method="kendall") 
nb_k <- expression("Kendall's"~tau == 0.964)

sp_nb <- ggplot(rankD, aes(x=r_Bentham, y=r_Nash))
sp_nb <- sp_nb + geom_point() + theme_bw() + 
  xlab(expression(paste("Bentham-based rank (", alpha %->% 1, ')')))+
  ylab(expression(paste("Nash-based rank (", alpha %->% 0, ')')))+
  annotate("text", x=22, y=70, parse = T, label=as.character(nb_s), size=4) +
  annotate("text", x=19, y=65, parse = T, label=as.character(nb_k), size=4) +
  labs(title="(c)")+theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.title = element_text(size = rel(1)))+
  theme(axis.title.x = element_text(size=12),axis.title.y = element_text(size=12)) 

#2 rank correlation, Spearman & Kendall, Hill & Webster
cor.test(rankD$Hill, rankD$Webster, method="spearman")
hw_s <- expression("Spearman's"~rho == 0.969)
cor.test(rankD$Hill, rankD$Webster, method="kendall") 
hw_k <- expression("Kendall's"~tau == "0.880")

#check labels
sp_hw1 <- ggplot(rankD, aes(x=r_Webster, y=r_Hill))
sp_hw1 <- sp_hw1 + geom_point() + theme_bw() + xlab("Webster-based rank") + 
  ylab("Hill-based rank") + annotate("text", x=22, y=70, parse = T, label=as.character(hw_s), size=9) +
  annotate("text", x=19, y=65, parse = T, label=as.character(hw_k), size=9) + 
  geom_text(aes(label=rownames(rankD)), size=4)
sp_hw1

sp_hw2 <- ggplot(rankD, aes(x=r_Webster, y=r_Hill))
sp_hw2 <- sp_hw2 + geom_point() + theme_bw() + 
  xlab(expression(paste("Webster-based rank (", alpha,' = 2)'))) + 
  ylab(expression(paste("Hill-based rank (", alpha,' = -1)'))) + 
  annotate("text", x=22, y=70, parse = T, label=as.character(hw_s), size=4) +
  annotate("text", x=19, y=65, parse = T, label=as.character(hw_k), size=4) + 
  labs(title="(d)")+theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.title = element_text(size = rel(1)))+
  theme(axis.title.x = element_text(size=12),axis.title.y = element_text(size=12)) 

#3 rank correlation, Spearman & Kendall, Adams & Jefferson
cor.test(rankD$Adams, rankD$Jefferson, method="spearman")
aj_s <- expression("Spearman's"~rho == 0.865)
cor.test(rankD$Adams, rankD$Jefferson, method="kendall") 
aj_k <- expression("Kendall's"~tau == 0.702)

#check labels
sp_aj1.1 <- ggplot(rankD, aes(x=r_Jefferson, y=r_Adams))
sp_aj1.1 <- sp_aj1.1 + geom_point() + theme_bw() + xlab("Jefferson-based rank") + 
  ylab("Adams-based rank") + annotate("text", x=22, y=70, parse = T, label=as.character(aj_s), size=9) +
  annotate("text", x=19, y=65, parse = T, label=as.character(aj_k), size=9) + 
  geom_text(aes(label=rownames(rankD)), size=4)
sp_aj1.1

sp_aj1.2 <- ggplot(rankD, aes(x=r_Jefferson, y=r_Adams))
sp_aj1.2 <- sp_aj1.2 + geom_point() + theme_bw() + 
  xlab(expression(paste("Jefferson-based rank (", alpha,' = 10)'))) + 
  ylab(expression(paste("Adams-based rank (", alpha,' = -9)'))) + 
  annotate("text", x=22, y=70, parse = T, label=as.character(aj_s), size=4) +
  annotate("text", x=19, y=65, parse = T, label=as.character(aj_k), size=4) + 
  annotate("text", x=19, y=50, label="Greece", size=3) +
  annotate("text", x=24, y=42, label="Denmark", size=3) +
  annotate("text", x=32, y=50, label="Jamaica", size=3) +
  annotate("text", x=38, y=9, label="Switzerland", size=3) +
  annotate("text", x=57, y=20, label="Paraguay", size=3) +
  annotate("text", x=65, y=32, label="Mali", size=3) +
  annotate("text", x=71, y=49, label="Russia", size=3) +
  labs(title="(e)")+theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.title = element_text(size = rel(1)))+
  theme(axis.title.x = element_text(size=12),axis.title.y = element_text(size=12)) 



# highlinging the 1996 Italian election in Samuels & Snyder because of their errors 

sp_aj2.1 <- ggplot(rankD, aes(x=r_Jefferson, y=r_Adams))
sp_aj2.1 <- sp_aj2.1 + geom_segment(aes(x = 4, y = 10, xend = 4, yend = 35)) + 
  geom_point() + theme_bw() + xlab(expression(paste("Jefferson-based rank (", alpha,' = 10)'))) + 
  ylab(expression(paste("Adams-based rank (", alpha,' = -9)'))) +
  annotate("text", x=10, y=40, label="Italy", size=2.7) +
  annotate("text", x=10, y=37.5, label="(Samuels & Snyder)", size=2.6) +
  annotate("point", x=17.5, y=20.33333, pch=17, colour= "red", cex=2) + 
  annotate("text", x=10.5, y=21, label="Italy1996", size=2.7) +
  annotate("point", x=18.5, y=20.66666, pch=17, colour= "red", cex=2) + 
  annotate("text", x=15.7, y=23.8, label="Italy2001", size=2.7) +
  annotate("point", x=20.33333, y=69.33333, pch=17, colour= "red", cex=2) + 
  annotate("text", x=18, y=67, label="Italy2006", size=2.7) +
  annotate("point", x=20.66666, y=69.66666, pch=17, colour= "red", cex=2) + 
  annotate("text", x=20, y=72, label="Italy2008", size=2.7) +
  annotate("point", x=69.75, y=26.5, pch=15, colour= "green", cex=2) +
  annotate("text", x=67, y=29, label="Taiwan1998", size=2.7) +
  annotate("point", x=69.5, y=23.66666, pch=15, colour= "green", cex=2) +
  annotate("text", x=60, y=25.3, label="Taiwan2001", size=2.7) + 
  annotate("point", x=69.25, y=23.33333, pch=15, colour= "green", cex=2) +
  annotate("text", x=65, y=20, label="Taiwan2004", size=2.7) +
  annotate("point", x=71.5, y=51.5, pch=15, colour= "green", cex=2) +
  annotate("text", x=66, y=49, label="Taiwan2008", size=2.7) +
  labs(title="(f)")+theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.title = element_text(size = rel(1)))+
  theme(axis.title.x = element_text(size=12),axis.title.y = element_text(size=12)) 

# combining graphs 0 to 4 (graph0:cma_mn2, graph1:sp_nb, graph2:sp_hw2, graph3:sp_aj1.2, graph4:sp_aj2.1)
# library(grid)
pdf("Figure3.pdf", width=7, height=11)

grid.arrange(cma_mn2, cma_mb2, sp_nb, sp_hw2, sp_aj1.2, sp_aj2.1, ncol=2)

dev.off()


################
### Figure 4  ##
################

# barplot
# This R script is for stacked bar chart.
#Italy
italy <- read.csv("a_0.csv", header=TRUE)
italy <- italy[,-6:-9]
summary(italy)
class(italy)
View(italy)
i <- cbind(
  ITA1996 = c(0.000000000, 0.000597263, 0.002509388),
  ITA2001 = c(0.000000000, 0.000494242, 0.002909469),
  ITA2006 = c(0.022336517, 0.002521247, 0.000000000),
  ITA2008 = c(0.026884233, 0.002129631, 0.000000000)
)

i


# labels
labels <- c("Abroad/Indigenous", "Apportionment", "Districting")

#Taiwan
taiwan <- read.csv("a_0.csv", header=TRUE)
taiwan <- taiwan[,c(1,6:9)]
summary(taiwan)
class(taiwan)
View(taiwan)
t <- cbind(
  TWN1998 = c(0.0124651890, 0.0078665720, 0.0000985071),
  TWN2001 = c(0.011938605, 0.007141397, 0.000107738),
  TWN2004 = c(0.011044676, 0.006732506, 0.000380146),
  TWN2008 = c(0.032624512, 0.030155742, 0.002028755)
)

t

# Italy & Taiwan
pdf("Figure4.pdf")
par(oma = c(2, 2, 0, 0)) #margin
par(mfrow = c(2,2))      # 2 by 2 matrix
par(mar = c(3, 3, 2, 1)) # the level of margin for x lines [lower, left, upper, right]
barplot(i, legend.text = labels, 
        args.legend = list(x = 3, y = max(colSums(t))),
        ylim=c(0,0.07), cex.names=0.8   
)
mtext(text=expression(paste(alpha,"-divergence (", alpha %->% 0, ')')),side = 2, line = 3)
barplot(t, legend.text = NULL,
        ylim=c(0,0.07), cex.names=0.8
)
barplot(
  sweep(i, 2, 100 / colSums(i), "*"), cex.names=0.8
)
mtext(text="%",side = 2, line = 3)
barplot(
  sweep(t, 2, 100 / colSums(t), "*"), cex.names=0.8
)

mtext("Election",
      outer = TRUE,      
      side = 1,         # lower side 
      cex = 1,         
      line = 0
) 
dev.off()

# check
# Italy & Taiwan
png(filename = "Figure4Check.png", width = 600, height = 600)
par(oma = c(2, 0, 2, 0))
par(mfrow = c(2,2))
par(mar = c(3, 3, 2, 1)) 
iprop = prop.table(i, margin = 2) #proportion
tprop = prop.table(t, margin = 2) #proportion
barplot(i, legend.text = labels, 
        args.legend = list(x = 2.4, y = max(colSums(t))),
        ylim=c(0,0.07)        
)
barplot(t, legend.text = NULL,
        ylim=c(0,0.07)
)
barplot(
  iprop
)

barplot(
  tprop
)
mtext("Real",
      outer = TRUE,      
      side = 3,          
      cex = 1.5,         
      line = 0
) 
mtext("%",
      outer = TRUE,      
      side = 1,          
      cex = 1.8,         
      line = 0
) 
dev.off()


################
### Figure 5  ##
################
# lineplot
# Italy (2008), graph 5

it08 <- read.csv("Italy2008.csv", header=TRUE)
colnames(it08)[2] <- "Population"
colnames(it08)[3] <- "Seats"
it08m <- melt(it08,
              id="district.name",
              measure=c(
                "Population",
                "Seats"
              ))
colnames(it08m)[1] <- "District"
colnames(it08m)[2] <- "Quotient"
it08m_l <- ggplot(data=it08m, aes(x=District, y=value, group = Quotient, linetype=Quotient))
it08m_l <- it08m_l + theme_bw() + geom_line() + xlab("Districts") + 
  ylab("Quotient")  +  
  labs(title="Italy (2008)", size=25) +
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0.3, size=16),
        axis.text.y=element_text(size=20))+
  theme(axis.title.x = element_text(size=18), axis.title.y = element_text(size=18)) +
  scale_x_discrete(limits=c("EMILIA-ROMAGNA", "PUGLIA", "LOMBARDIA2", "LAZIO1", "LOMBARDIA1", "TOSCANA", "CAMPANIA1", "VENETO1", "CAMPANIA2", "SICILIA2", "SICILIA1", "PIEMONTE1", "PIEMONTE2", "EUROPA", "CALABRIA", "VENETO2", "SARDEGNA", "LIGURIA", "LAZIO2", "MARCHE", "LOMBARDIA3", "ABRUZZO", "FRIULI-VENEZIA GIULIA", "AMERICA MERIDIONALE", "TRENTINO-ALTO ADIGE", "UMBRIA", "BASILICATA", "AMERICA SETTENTRIONALE E CENTRALE", "MOLISE", "AFRICA ASIA OCEANIA ANTARTIDE", "VALLE D'AOSTA"
  )
  )+
  geom_vline(xintercept = 14, colour="red") + 
  geom_vline(xintercept = 24, colour="red") +
  geom_vline(xintercept = 28, colour="red") +
  geom_vline(xintercept = 30, colour="red") + 
  theme(legend.position="none")+
  theme(plot.title = element_text(size = rel(2)))
print(it08m_l)

# Taiwan (2004), graph 6
ta04 <- read.csv("Taiwan2004.csv", header=TRUE) # Aborigine --> Indigineous
colnames(ta04)[2] <- "Population"
colnames(ta04)[3] <- "Seats"
ta04m <- melt(ta04,
              id="district.name",
              measure=c(
                "Population",
                "Seats"
              ))
colnames(ta04m)[1] <- "District"
colnames(ta04m)[2] <- "Quotient"
ta04m_l<- ggplot(data=ta04m, aes(x=District, y=value, group = Quotient, linetype=Quotient))
ta04m_l <- ta04m_l + theme_bw() + geom_line() + xlab("Districts") + 
  ylab("Quotient")  +  
  labs(title="Taiwan (2004)", size=25) + 
  theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0.3, size=20),
        axis.text.y=element_text(size=20))+
  theme(axis.title.x = element_text(size=18), axis.title.y = element_text(size=18))+
  scale_x_discrete(limits=c("Taoyuan County", "Taichung County", "Taipei City 1", "Taipei County 2", "Taipei City 2", "Changhua County", "Taipei County 3", "Kaohsiung County", "Tainan County", "Taipei County 1", "Taichung City", "Kaohsiung City 1", "Pingtung County", "Yunlin County", "Tainan City", "Kaohsiung City 2", "Chiayi County", "Miaoli County", "Nantou County", "Yilan County", "Hsinchu County", "Keelung City", "Hsinchu City", "Hualien County", "Chiayi City", "Indigenous (Mountain)", "Indigenous (Plain)", "Taitung County", "Penghu County", "Kinmen County", "Lienchiang County"
  )) + 
  geom_vline(xintercept = 26, colour="red") + 
  geom_vline(xintercept = 27, colour="red") + 
  theme(legend.position="none")+
  theme(plot.title = element_text(size = rel(2))) +
  theme(plot.title = element_text(hjust = 0.5))
print(ta04m_l)

# Taiwan (2008), graph 7
ta08 <- read.csv("Taiwan2008.csv", header=TRUE) # Aborigine --> Indigineous
colnames(ta08)[2] <- "Population"
colnames(ta08)[3] <- "Seats"
ta08m <- melt(ta08,
              id="district.name",
              measure=c(
                "Population",
                "Seats"
              ))
colnames(ta08m)[1] <- "District"
colnames(ta08m)[2] <- "Quotient"
ta08m_l<- ggplot(data=ta08m, aes(x=District, y=value, group = Quotient, linetype=Quotient))
ta08m_l <- ta08m_l + theme_bw() + geom_line() + xlab("Districts") + 
  ylab("Quotient")  +  
  labs(title="Taiwan (2008)", size=25) +
  theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0.3, size=16),
        axis.text.y=element_text(size=15))+
  theme(axis.title.x = element_text(size=18), axis.title.y = element_text(size=18)) +
  scale_x_discrete(limits=c("Hsinchu County", "Yilan County", "Keelung City", "Tainan City 1", "Tainan County 3", "Hsinchu City", "Tainan County 2", "Taichung City 2", "Tainan City 2", "Taipei City 4", "Yunlin County 2", "Yunlin County 1", "Tainan County 1", "Taipei City 3", "Changhua County 3", "Taipei County 8", "Taipei City 1", "Kaohsiung City 1", "Taipei County 11", "Taipei County 4", "Taichung County 2", "Taipei City 2", "Kaohsiung County 4", "Changhua County 4", "Taichung County 3", "Taipei County 1", "Taipei County 3", "Kaohsiung City 4", "Kaohsiung County 2", "Taipei City 6", "Taipei City 5", "Taipei City 7", "Taichung City 1", "Kaohsiung City 5", "Taipei County 2", "Taichung City 3", "Taipei County 9", "Taipei County 10", "Taoyuan County 2", "Taipei City 8", "Changhua County 1", "Taoyuan County 1", "Taipei County 12", "Taoyuan County 4", "Changhua County 2", "Taoyuan County 3", "Pingtung County 1", "Kaohsiung County 1", "Taoyuan County 5", "Chiayi County 2", "Kaohsiung County 3", "Taoyuan County 6", "Taichung County 5", "Miaoli County 2", "Kaohsiung City 2", "Pingtung County 2", "Taipei County 5", "Taipei County 6", "Taipei County 7", "Pingtung County 3", "Chiayi County 1", "Nantou County 2", "Chiayi City", "Taichung County 4", "Taichung County 1", "Kaohsiung City 3", "Hualien County", "Miaoli County 1", "Nantou County 1", "Indigenous (Mountain)", "Indigenous (Plain)", "Taitung County", "Penghu County", "Kinmen County", "Lienchiang County"
  )) + 
  geom_vline(xintercept = 70, colour="red") + 
  geom_vline(xintercept = 71, colour="red") + 
  theme(legend.position="none")+
  theme(plot.title = element_text(size = rel(2)))+
  theme(plot.title = element_text(hjust = 0.5))
print(ta08m_l)


# combining graphs 5 to 7
f5 <- grid.arrange(it08m_l, arrangeGrob(ta04m_l, ta08m_l, nrow=2), heights=c(2/5, 3/5), ncol=1)
ggsave(filename="Figure5.pdf",plot=f5, units="in", width=20,height=20)
